load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            User Config                   ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

configFile = "StationLists/YLJ_precip_stnlist.txt"
lines = readAsciiTable(configFile,1,"string",1)
delim = ","

stnlist   = ndtooned(str_get_field(lines,1,delim))
stnnum    = ndtooned(str_get_field(lines,2,delim))
stnlat    = ndtooned(tofloat(str_get_field(lines,3,delim)))
stnlon    = ndtooned(tofloat(str_get_field(lines,4,delim)))
stndomain = ndtooned(str_get_field(lines,5,delim))
nstns     = dimsizes(stnlist)

precip_dt = 1 ; output 1h precip

if (.not. isvar("input_file")) then
    input_file = "/public/users/wrf/wrfout/2012081808/wrfout_d04_2012-08-18_00:00:00"
end if
if (.not. isvar("output_prefix ")) then
    output_prefix = "~/test/"
end if
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;           Open NetCDF File               ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
a = addfile(input_file + ".nc", "r")
filename = systemfunc("basename " + input_file)

domain = substring(filename, 7, 9)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;              Process Time                ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times = a->Times
hours = WRF_Times2Udunits_c(times, 0)
wrf_dt = hours(1) - hours(0)
; UTC -> CST
hours = hours + 8 ; 
timeStr_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")
ntimes = dimsizes(hours)

dt = tointeger(precip_dt / wrf_dt + 0.1)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            Read Precipitation            ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
rain_exp = wrf_user_getvar(a,"RAINNC",-1)
rain_con = wrf_user_getvar(a,"RAINC",-1)
rain_tot = rain_exp + rain_con

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;          Loop for each station           ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Setting Output Filename
output_filename = output_prefix + "/precip_output_YLJ_" +  "_" + ut_string(hours(0), "%Y-%N-%D_%H") + ".txt"
print("Creating: " + output_filename)

do istn = 0, nstns - 1, 1
    if (.not. isStrSubset(filename, stndomain(istn))) then
        continue
    end if

    ;;;;;;;;;;;;;;;;;; latlon to XY ;;;;;;;;;;;;;;;;;;
    res = True
    res@returnInt = False   ; return real values
    loc = wrf_user_ll_to_ij(a, stnlon(istn), stnlat(istn), res) - 1

    syc0 = floattointeger(loc(0))
    sxc0 = floattointeger(loc(1))

    syc1 = floattointeger(loc(0)) + 1
    sxc1 = floattointeger(loc(1)) + 1

    yy = (/syc0, syc1/)
    xx = (/sxc0, sxc1/)

    print ("X/Y location is: " + loc)

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; crop station area from total rain data
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    rain_tot_crop = rain_tot(:, syc0:syc1, sxc0:sxc1)
    
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; calculate interpolated stational rain
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    rain_tot_station = ndtooned(linint2(xx, yy, rain_tot_crop, False, loc(1), loc(0), 0))

    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; output precipitation tendency for dt hours
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    line_array = new( tointeger(ntimes/dt) + 1, "string" )
    l_cnt = 0
    ; write header line
    line_array(l_cnt) = "Station Number, Station Name, Start Time, End Time, Precipitation"
    l_cnt = l_cnt + 1
    do it = dt, ntimes-1, dt   ; TIME LOOP
        rain_tend = rain_tot_station(it) - rain_tot_station(it - dt)
        line = stnnum(istn) +", "+ stnlist(istn) +", "+ timeStr_cst(it - dt) + ", " + timeStr_cst(it) + ", " + sprintf("%5.2f", rain_tend)
        line_array(l_cnt) = line
        l_cnt = l_cnt + 1
    end do
    asciiwrite(output_filename, line_array)
end do
